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Abstract 

In the paper we describe a superexponentially convergent numerical-analytical method for solving 
the eigenvalue problem for the some class of singular differential operators with boundary conditions. 

^^ The method (FD-method) was firstly proposed by V. L. Makarov and successfully combines the 

benefits of using the coefficient approximation methods (CAM) and the homotopy approach. The 
sufficient conditions which provides convergence of the proposed method are stated and rigorously 
substantiated. The algorithm for the software implementation of the proposed method is described 
too. A lot of numerical examples are included in the paper. The examples confirm the theoretical 
conclusions. We also have made the comparison between the results obtained by FD-method and 

'^ ' results obtained by the powerful software package for solving Sturm-Liouville problems — SLEIGN2. 

^ MSC 2010: 65L15, 65L20, 33D15, 68W99 

g 1 Introduction 

,__! There is a great scope of the numerical methods for solving eigenvalue problems for 

^ differential operators of the second order (see pj). The known numerical techniques for 

0\ eigenvalue problems can be divided into two groups: methods based on the direct approx- 

00 

l/~j imation of the solution of differential equation and methods based on approximation of 

I^ the coefficients of differential equation (see [T]). 

,-H The methods of the first group such as finite difference, finite elements and spectral are 

extensively treated through both theoretical investigations and well developed soft. The 
main idea of this approach is replacement of the eigenfunctions by piecewise polynomial 
functions, which results in the approximation of a differential equation by a system of 
linear algebraic equations. Because of the different nature of the original differential 
operator and approximating algebraic operator the numerical solution comes close to the 
exact one only for the low-indexed eigenvalues and corresponding eigenfunctions, that is 
defined by the size of the mesh step. Hence, such numerical techniques are effective to 
find the low-indexed eigenvalues, but not effective to find the high-indexed eigenvalues [T] . 
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The methods of the second group are based on the general idea of the coefficient ap- 
proximation methods (hereinafter referred to as CAM), also known as the Pruess methods, 
— to replace the coefficient functions of the problem by simpler approximation and then 
to solve the approximating problem (see J9]). Firstly this idea was proposed by Kryloff 
and Bogolioubov (1926) in J4] for the case of piecewise constant approximation. S. Pruess 
provided a thorough convergence and error analysis of such methods when piecewise poly- 
nomial approximation is used. Except for Gordon, who uses linear functions, the earlier 
references all confine themselves to piecewise constant approximations, and this is the 
most practical choice: otherwise the approximating problem may be no easier to solve 
numerically than is the original. It is well known (see pi) that if the piecewise constant 
approximation with the maximum step size h is used, the error of the CAM is estimated 
by 0{h). Hence, when we using the CAM, a crucial problem arise: to halve the error we 
have to double the computational time. This problem in many cases can be solved by 
using the functional- discrete method (hereinafter referred to as the FD-method) which is 
the essential generalization of the CAM for solving Sturm-Liouville problems of different 
types. The algorithm of FD-method was firstly described and justified for linear Sturm- 
Liouville problem in [6] . Generally speaking, it consists of two main steps: the first one is 
to find the initial approximation by applying the CAM and the second one is to calculate 
a sufficient number of corrections to achieve needed accuracy. The second step is substan- 
tiated by the homotopy method (also known as perturbation method, continuation method 
or successive loading method, see, for example, [2]). From that point of view, FD-method 
is the synthesis of the CAM and homotopy method. 

In the present paper we consider the following singular Sturm-Liouville problem 

q{x)u{x) = Xu{x), xG(— 1,1), (1.1) 



dx 



(1 _ a;2^) M^ 



dx 



lim (i_a;2)^^=o. (1.2) 

x^±i ^ ' dx 

Since the problem (1.1), (1.2) has singularities at the points ±1, we can't directly apply 
the algorithm of the FD-method, proposed in (gI. However, in |5] the modified FD- 
method technique for the problems with such singularities have been presented. Let us 
state briefly the algorithm of the FD-method described in [5]. According to the FD- 



method we approximate the eigenvalues and eigenfunctions of the problem (1.1), (1.2) by 
the following truncated series 



A„=5^A(/), u^{x) = Y,^^\^)^ (1-3) 

j=0 i=0 

where the unknown summands Xt and Un \x), for i = 1,2,..., rn^H are the solutions to 



(^The positive integer number m is called the rank of FD-method. 



the following recurrence problems: 



d 
dx 



(1-x^) 



dUn {x) 



dx 



i-1 



+ A»4^-)(x) = -5^A(f-)«i^)(x) + g(x)4^--i)(x) = Fi^-)(x), (1.4) 



s=0 



X G (— 1, 1) with the boundary conditions 



lim 



[l-X^) 



.U) 



du\V{x) 
dx 



1.5) 



and the orthogonality conditions 



u^^\x)u^^\x)dx = 0, j = 1,2, . . . ,m. 



;i.6) 



-1 



The initial data A™ , Un {x) which is needed to begin the recurrence process can be 
found as the solution to the following eigenvalue problem, called the basic problem: 

dUn (x) 



d 

dx 



(1 - x^) 



dx 



lim 



+ Ai%f(x) = 0, XG(-1,1) 

dUn\x) 



(1 - X^) 



dx 



0. 



The solutions to problem (1.7), (1.8) can be expressed in the following form 



AW=n(n + l), «f(x) = y^ 



2n + l 



-* n \X ) , 



;i-7) 



(1- 



;i.9) 



where Pn{x) denotes the Legendre polynomial of the first kind of order n G N (see [3j). 
The following statement about the convergence of the FD-method, described above, was 
proved in JSJ. 

Theorem 1.1 Let the function q{x) belongs to the Q^ [—1, 1] , — the space of piecewise 
continuous functions defined on [—1, 1] . If there exists a nonnegative integer no, such that 

2|k(2;)|loo,hi,i] 



* 



no 



no 



< 1, no > 1, go = 2||g(x)||^_j_^_^] < 1 



1.10) 



then the FD-method (1.4) -(1.9) has a superexponential convergence rate for all numbers 
n > Ho : 



Un{x)- Un [X) < 



n = no, no + 1, ... , where 



(Qu 



,m+l 



1 -gn 



(2m -1)!! 

«m+l — 2- TTT < 



■«m+l, 



1 



m 
■^n An 






1.11) 



(2m + 2)!! ~ {m + l)^/7^m' 



lk(^)lloo.[-i,i] = max |g(x)|. 

xg — 1,1 



Theorem |1.1| was proved in [5] for the case when q{x) = q{—x) but, using the same 
technique, it can be proved for the case of an arbitrary function from Q^ [—1, 1] without 
considerable difficulties. 



If we try to apply the FD-method (1.4) -(1.9) to approximate eigensolution A„,u„(a; 



with n = 0, 1, . . . , no — 1 we might find it divergent. The reason of that is the lack of ap- 



proximation for q(x) in equation (1.7) of the basic problem. To overcome such difficulties 



we have to use the general scheme of the FD-method which uses the approximation of the 
function q{x) on [—1, 1] by the piecewise constant function q{x). 



2 The general algorithm of the FD-method for solving the 
Sturm-Liouville problem with Legendre differential operator: 
theoretical justification 

To make the general algorithm of the FD-method more easy for understanding, let us 
consider the following general problem 

d_ I. ^2. du{x,t) 
dx dx 



+ [\{t) - q{x) - t {q{x) - q{x))] m(x, t) = 0, x G (-1, 1) , t G [0, 1] 

(2.1) 
(2.2) 



lim (l_a;2)^!^ = 0, Vt G [0, 1] 



dx 
To define the function g(x), we have to fix some mesh on the interval [—1,1] : 

uj = | — 1 = xn < a;i < . . . < xn^i < a^Af = 1| , h = max (xi — Xi-i) . 

^ ' i=l,2,...,N 



(2.3) 



We also require that all points of discontinuity of the function q{x) on the interval [—1, 1] 
are contained in the set u. 

There are many possibilities to define q{x), for example, 

'1 



q{x) = q{x,uj) 



q [ - {Xi_i + Xi 



X G IXj_i, Xi 



1,2,. ..,iV, 



or 



q{x) = q{x, w) = 2 ("Jl^i-O + Qi^i)) ^ ^ ^ [xi^i,Xi) , 



2 = 1,2,. ..,iV, (2.4) 

but, in any case, the piecewise constant function q{x, u) must have the following property 

0. 



lim I max \q(x,uj) — q(x) 
h^o V xe[-i,i] 



It is easy to see that if we set t = 1, then problem (2.1), (2.2) would transform to the 



problem (|l.l|), (|1.2|). Assume that the solution to the problem (2.1), (2.2) can be expressed 

oo oo 



in the series form 



[X 



(2.5) 



j=0 



j=0 



and 

dun{x,t) 
dx 



E'^'f^. "^-f'''^. V..(-U).W.|0.1]. 



i=o 



j=o 



dx'^ 



(2.6) 

The unknown summands of the series (2.5 ) can be found as the solutions to the recurrence 
problems (j = 0, 1, . . . , m): 



d 
dx 



(1-x^) 



dui\x) 
dx 



+ [AW-g(x)]t.(f)(a:) = 

x), xG(-1, 1), 



(2.7) 



with the boundary conditions 



lim 



matching conditions 

and additional requirement 



(1-x^) 



dUn {x) 

dx 



dUn {x) 

dx 



0, 



(2. 



0,i = 1,2,...,A^-1. 



(2.9) 



u:fl' {x)u\' {x)dx = 6oj 



-1 



1, ifj = 0, 
0, ifjVO. 



(2.10) 



Here square brackets denotes the jump of the function at the point x = Xi. 

Among the problems (2.7) - (|2.9) it is worth to emphasize the first one (j = 0): 



dx 



:i-x 



lim 

X-¥±l 



2.dUn (x) 



dx 



{1-X-) 



+ [Af -g(x)]tzf(x) = 0, xG(-l,l), 

dUn {x) 



dx 



{u^^\x)) dx = l, 



Wna:)]_^ = o, 



dUn (x) 

dx 



0. 



(2.11) 



(2.12) 



(2.13) 



easy to make sure that the differential operator 



As it was mentioned above, the problem (2.11) - (2.13) is called the basic problem. It is 

operato 

(l-x^)^ -g(x)(-) 



L ■] 



d 
dx 



dx 



(2.14) 



is self-adjoint in the Sobolev space 



ly'"' (-1, 1) = h{x) e W'' (-1, 1) : ^Um (1 - x') fix) = o| . 



.(0) 



This fact imphes that the eigenfunctions Un (x), n = 0, 1, . . . of the problem (2.11) - (2.13) 



■2,1 



form the complete orthogonal system in W ' (—1, 1) and the corresponding eigenvalues 
An , (Aj- < Xj , when i < j, i, j = 0,1,2 . . .) are simple. 

Suppose that the basic problem is solved and a pare Xn ,Un (x) denotes its arbitrary 
eigensolution. 



It is well known that the problems (2.7) - (2.9) for j = 1,2, ... ,711 can be solved if and 
only if the following equality holds 

1 



u^^\x)Fi^\x)dx = 0. 



(2.15) 



Equality (2.15) together with (2.7) yield us the formula for computing A 



n 



^i''^ = -Y.^t'^ I ''n\^Wn\^)dx + I {q{x) - q{x))uli-'\x)u^^\x)dx. (2.16) 



Taking into account (2.10) we can simplify formula (2.16) in the following way 

1 



A(f) = / {q{x)-q{x))u^r'\^Wn\x)dx. 



-1 



,(J) 



The function Un (x) can be represented as a Fourier series 



(2.17) 



ul^){x) = 5^uf (x) / ul^\x)uf\x)dx 



s=0 



(2.18) 



To find the Fourier coefficients let us multiply both sides of the equation (2.7) by ui (x) 
and then integrate them on the interval [—1, 1] : 



1 



X] 



d 
dx 



{1-x^) 



dui\x) 
dx 



dx+ u'f\x) [Af - q{x)] uli\x)dx = (2.19) 



u^^\x)Fi^\x)dx. 



Taking into account boundary conditions (2.8) and using the integration by parts (two 



times), from (2.19) we obtain 



ul^\x) 



dx 



(l-xO 



dus (x) 
dx 



q[x)ui^\x)u^^\x) I dx + A(f) / ul^\x)u^^\x)dx = 

(2.20) 



j uf\x)Fi'\x)dx. 



And (2.20) immediately lead us to the equality 



^i3)(^.\„M 



u^i\x)u^^\x)dx 



-1 



/ Fi^'\x)ur{x)dx 



x(o) x(0) 
An — As 



(2.21) 



Using (2.21) we can express the formula (2.19) in the following form 

1 



u'^JHx) 



^ j Fk'\x)uf{x)dx 
2-^ \(0) _ x(o) -p 

p=0 "^n Ap 



u^'\x). 



(2.22) 



Formulas (2.17) and (2.22) yields the estimations 



lAi^^l < \\q{x) -q{x) 



loo, [-1,1] IP" (^) 



(2.23) 



and 



\uii\x)\\ < Mn\\F}^\x)\\ = M„ 



•j-i 



i-i 

s=l 



(2.24) 



<MAY1 \>^n-'^\ \H\^)\\ + h(^) - ^(^)lloo,[-l,l] H~'\^)\\ < 



. s=l 



i-1 



< M„ \\q{x) - g(x)lL^j_,^i] <^ Yl hn'' 

I s=0 

respectively, where 

M. = max(rAf-Ail)",rAi°|,-Af 






n = l,2,..., (2.25) 



Mo = (Af ) - A, 



(0) _ x(0)\ ^ ^ 



and ||-|| denotes the common L2-norm, ||/(a;)|| = \ J f'^{x)dx. 



using estimate (2.23), we obtain 



Multiplying both sides of the inequality (2.24) by ( M„ \\(l{'^) ~ ?(^)lloo f-i ii 



ui' ix) 



i-1 



Ur, 



\x) 



Mn ix] 



M„ ||g(x) - g(x)|L^f_i_i])' s=o [Mn \\q{x) - q{x)\\^^^_^^^^ 
Using the notation 



j-i' 



and 



(2.26) 



Vj = \\uii\x) 



(m„ 



\q[x) — q[x) 



loo, [-1,1] 



ui^Hx)\\r~ 



(2.27) 



we can represent inequality (2.26) in the following form 

i-i 



^ ^"^j-s-iVsJ = 1,2,..., vo = \\ul^\x)\\ = 1. 



(2.28) 



s=0 



It is easy to see that the sequence {Vi} , defined by the following recurrence formula. 



i-i 



^J=5Z^^-i-^^' ^0 = 1, J = 0,1,... 



s=0 



satisfies inequalities 

Vj > Vj, j = 0, 1, . . . . 

Let us consider the generating function f{z) defined by the formula 



fi^) = J2^s 



z\ zeK. 



s=0 



From (2.29) we obtain the equation with respect to f{z) 

f{z) = zf{z) + l, ze{-p,p), 



(2.29) 



(2.30) 



(2.31) 



(2.32) 



where p denotes the convergence radius for the power series (2.31). The solution to 



equation (2.32), which satisfies the condition /(O) = 1 is 



fiz) = (2z)-i (1 - VT^^z) = ( 2^ + f; ^^'iJ^'h ^zY 1 i2z)-\ z e 



s=2 



(2.) 



1 1 

"4'4 



and 



' (2j+2)!! 



dof 



4^ = «.4^ J = 1,2,..., (-1)!! = 1. 



(2.27), (2.30) together with (2.34) lead us to the estimation 

" ^'\x)\\ = (Mu Hx) - qix)\\^^^_,y Vj < {AfnY aj = rlaj. 



Wl: 



(2.33) 
(2.34) 

(2.35) 
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And inequality (2.35) yields the following error estimations of the FD-method's conver- 
gence rate 



Unix)- Un (x) 



< < 



^■m+l 



1 — r, 



-«m+l, 



■^n Ar 



< < 



\\q{x) — q{x) 
\\q{x) -q{x) 



E ^= 

, j=m+l (j + 1)V7U 



if r„ < 1 



, if rn 



(2.36) 



-ar 



1 

j~m {] + 1)0U 



(2.37) 



OD 
loo, [-1,1] 1^ 

j= 
The infinite sum in the right sides of inequalities (|2.36[) and (2.37) can be easily estimated 



if rn < 1, 
, if r„ = 1. 



E 



< 



,_™(j + l)v^ J ix + l)y/xTT 

■?-'" m-1 



dx 



2 / 1 

^= arctan —^^=^— . . 
vr \ V m + 1 y 



(2.38) 



We have proved the following theorem. 
Theorem 2.1 Let q{x) G Q^[—l, 1] anc? i/ie constant fn = Mn \\q{x) — ^(x)||^ r i ii sat- 



isfies the inequality Vn < j. Then the FD-method (2.5), (2.7) - (2.9) converges faster 



than the geometric series with denominator r„ (superexponentially) . The estimations of 



its convergence rate are given by formulas (2.36), (2.37). 



Before passing to the next section, let us investigate the question about the asymptot- 



ical behavior of parameter M„ (2.25) as n tends to +oo. It is easy to see that the basic 



problem (2.11[) - (2.13) is equivalent to the following auxiliary eigenvalue problem 



d_ 
dx 



:i-x' 



dv{x,T) 
dx 



lim 



(1 - X^) 



+ [fi{r) - rq{x)] v{x, r) = 0, x G (-1, 1), r G [0, 1], (2.39) 

1 

dv(x,T) 



[v{x,T)i 



dx 

dv{x,T) 
dx 



0, I {v{x,T)fdx = l,WTe[0,l], 



0, VrG[0,l], z = l,2,...,iV-l, 



(2.40) 



(2.41) 



when r = 1. On the other hand, as it follows from the proof of Theorem |1.1| (see [5]), for 
sufficiently large number n, namely 



^>2|k(a:)IL,[_i,i]>2||g(x)||^^[_,^y, 



(2.42) 



theeigensolution/i„(r), f„(a;, r) to the problem (2.39) - (2.41) exists and can be developed 



to the power series with respect to r on the interval [0, 1], for all x G (— 1, 1), furthermore 

/i„(0)=n(ri + l),/i„(l) = A(f). (2.43) 



Hence, from equality (2.39) it is easy to obtain 
d 



dx 



l-x' 



d^Vn{x,T) 

dxdr 



d 

+ [/^n(i-) - rq{x)] -7^Vn{x, t) = [fx'^{r) - q{x)] Vn{x, r), (2.44) 



X G (—1,1), r G [0,1]. After that, applying to the both sides of equality (2.44) the 
1 
integral operator J Vn{x,T){-)dx, we obtain 



/^n(^) = / K(a;,7-)) q{x)dx, rG [0,1]. 



-1 



Taking into account (2.43) and equality (2.45), it is easy to obtain 



1 1 




A^ =n{n + l)+ / / {vn{x, T)f q{x)dxdT. 



-1 



Finally, equality (2.46) yields the following estimation 

Ai°|i-Af >2(n + l)-2||g(x)|L_[_,^,j>0 



(2.45) 



(2.46) 



(2.47) 



Hence, for the n sufficiently large (see (2.42)) we have the following estimation which 

describes the asymptotical behavior of the M„ : 

-1 



M„< (2n-2||g(x)|L^[_,^y 



(2.48) 



3 General algorithm of the FD-method: software implementa- 
tion. 



The solution to the basic problem (2.11) - (2.13) can be represented in the following 
form 



u^^>{x) = AP,^{x) + BiQ,^{x), xG [x,_i,Xi], A,B,eR, t = l,N, (3.1) 



where 



Ui = u. 



(Mf) 



i (-1 ± yi + 4(A--,)) 



(2 



qi = q{x), X G [xi_i,Xi) (3.2) 



and Qui{x) denotes the Legendre function of the second kind (see [3]). It is well known 
that (see, for example, (3)) 



(3.3) 



dPu{x) 2sin(7ri/) dPu{x) 
lim (1 — xn — ; = , lim (1 — x ) — ; = 0, 

x^-i dx IT 2.'->i dx 



dQ,{x) dQ,{x) 

iim (1 — X ) = cos (TTu) , hm (1 — x ) = 1. 

x^-i dx ^ ' ' x^i^ dx 



(^The sing "+" or "— " can be chosen optionally. 
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To satisfy boundary condition (2.12) we have to set 



2 sin (ttz/i) 
Ai hfiicos(7rz/i) = 0, Sat = 0. 



TT 



(3.4) 



The other constants A^, Bi have to be found from the matching conditions (2.13) as the 
solutions to the following recurrence systems of linear algebraic equations: 



= 2,3,. ..,N. 
Let us consider the function 



AiP^X^i^i) + BiQ^X^i_i), 
A,P;(x,_i) + B,g',^(a;,_i), 



(3.5) 



2 sin (vrz^i) 
$ (A) = Ai h Bi cos (vrz/i) 



TT 



where constants Ai,Bi are defined by the following recurrence formulas (see (3.5)) 

B,^i = — (AP.,_,(a;.-i) - C,P;_^(x,_i 



(5i_i = Py^^^{xi^i)Q'^^_^{xi-i) - P^^_^(xi_i)Q^,„,(xi_i), 

Ci = AiP^^Xxi^i) + BiQ^X^i^i), 

Di = AP'ixi-i) + B,Q'{xi.i), t = 2,3,...,N 



(3.6) 



with the initial data 



A 



N 



1, Bm = 0, 



(3.7) 



and the values u^ = z/j (A) are defined by formula (3.2 ). It is easy to see that the eigenvalues 



of the basic problem (2.11) - (2.13) coincide with the solutions of the equation 

$(A) = 0. 



(3.^ 



Since the eigenvalues of the basic problem are simple, they can be found with any precision 



as the solutions to equation (3.8) by the bisection method. 

Let An be an arbitrary solution to equation (3.8). Using the coefficients Ai,Bi,i 



AO) 



1,2, ... ,N, computed by formulas (3.6) with A = An we can compute corresponding 
eigenfunction uh, (x) by formula (3.1). 



,(i) 



To compute An we can use the formula 



Al^') 



1 -1 



(u^^'ix)) dx 



{q{x) — q{x)) u)i ^' {x)u'^' {x)dx . 



(3.9) 
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And for the computation of the functions Un (x), j = 1,2, . . . ,m it is convenient to use 
the formula 



v^^ix) = E,u'^\x) + / K{x,OFi^\Od^ 



(3.10) 



-1 -1 



where 



K{X,0 = Wn{x)u^^\0 - Wrr{0<K^) 



(0)/ 



Ej = - 



{<\x)ydx 



-1 



4°) (a;) K{x,OFi'\Od^\ dx 



-1 



-1 



and Wn{x) denotes a function, defined by the formula 



Wn{x) = AiPi^Xx) + BiQi,Xx), X G [xi-i,Xi] , i = 1,N, 



where Ai,Bi, i = 1,2,...,N can be computed recursively by formulas (3.6) with the 
initial data 

An = 0, Bn = 1. (3.11) 



Taking into account the properties of the Legendre functions (see [3]), it is easy to verify 
that 



dK{x,0 



dx 



(l-x^) 



2\-l 



:rG(-l,ll 



(3.12) 



^=x 



and K{x,^) is the Cauchy operator for the nonhomogeneous equations (2.7). 



In general case, the integrals in formulas (3.9) and (3.10) could not be calculated 



analytically. Moreover, the integrands in this formulas are unbounded at the points ±1. 
To compute this integrals it is convenient to use numerical methods, for example, the 
tank rule (see |8]) and Stenger's formula (see ^lOj): 

+00 



f{x)dx = f 



a + be'^ 



a) duj 



K 



risinc 2_^ J I ^ _L ^Ih 
l=-K ^ 



1 + e"^ / (e-^/2 + e^/2)' 
(6 -a) 



(3.13) 



+ e'' 



and 



2fc 



l=-K 



(p tflsinc/ -^ — L pj'i^sinc / '^\ 

ib-a) 



f(x)dx--h- V s^-^U /^ Q + ^g'^°" 

f[X) X ~ sine 2_^ k-l J y^j^ ^ih,,^^ J ^^_ih^^^j2 _^ ^lh,,„j2y 



(3 



(3.14) 



(^ The function f(x) is required to be sufficiently smooth on (a, b), see [lOl. 
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where 



Zk= , , _, . ,, , k = -K ...,K, 6^ 



(-1) 



1 



1 + e 



h ■ h 



+ 



sin (vrt) 



2 J 7Tt 





dt, k = -2K ...2K, 



and /lo 



277(4 

K 



Henceforth we require that the function q{x) is analytical on each interval (xj_i,Xj), 
i = 1,2, ... ,N. Before passing to the algorithm description, we need to introduce the 
following auxiliary notation 



Xj_i + Xie'^"^'^"^ 



l^i 



[•^i—l •^i) 



'*'^ l+g/i.^nci ' ^*'^ (^e-jhs,nc/2 _^f>jh,,„c/2'j'^' " 



-, -1 



{u^^\x)Ydx 



-1 



(3.15) 
t = l,2,. ..,N,j = -K,...,K. 

(k) (k) / 

One of the possible algorithms for computation of An and Un [k) for k = 1,2, ... ,r, 
i = 1,2, . . . ,N, j = —K . . . ,K is described below. 



Algorithm 1: The algorithm for computation of AJ^, u'}^{zij), d E l,r, i E 1,N, j E —K, K 



input : The arrays of values Un [zi^], Wn{zij), jj^ij, i G 1,N, j G —K,K, 






/ e -2K..2K, A°, and r — the depth of the FD-method. 



output: The arrays of values Un (zij), d ^ l,r, i ^ 1, N, j ^ —K, K, 



\^n\ del,r. 



begin 



for d := 1 to r do 



.Wi 



,(0) 



CoinputeLa.inhd8L{ Un (zij), Un (zij))// This procedure computes A 
see description below 



(d) 

n , 



M)l 



Ak) 



ConiputeF(Mn (^^jj), Mi)// This procedure computes the values 



Ad) 



Fn (zij) , see description below 



(0), 



M, 



ConiputeCorrectionForEigenfunction(M„ (2;ij), Wn{zij), Fn (zij)) 
// This procedure computes the values Un (zij), see description 



below 



id). 



.(0), 



ProvideOrthogonalityCondition( Un (zij), Un (zij))// This procedure 
recompute the values uh, (zij) to provide the orthogonality 

condition J Un {x)un {x)dx = 

-1 



end 
end 



(* About the optimal choice of the value for h^inc see lOl. 
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id). 



,(0)^ 



Procedure ConiputeLanibda( Un (zjj), u„ (zi.j)) 



,{d) 



(0), 



input : The arrays of values Un'{zij), Un (zij), i e l,iV, j <E —K,K 

output: The value A 

begin 



id) 



Xlf' :^ 



for i := 1 to A^ do 

for j := —K to K do 

Ad) ._ .(d) , ^ (0), 



.{d-l) 



id) 



A„ := A„ +Un'{zij)un {zi,j){qizij) ~ q{zij)} iiij // Compute A„ by formul 



(3.91, using formula (|3.13[) 



end 



end 

An '■= hginci^nMi // The An is computed 



end 





Procedure ProvideOrthogonalityCondition( Un (zij), Un \zij)) 


inpi 
out 




Lit : The arrays of values Un'[zi^j), Un^Zi^), i G 1, N, j E —K,K 


put: The array of values Un (zi^j), i E 1,N, j e —K, K, such that the integral 


/ Un{x)un{x), computed by formula (3.13), is equal to zero 
-1 


begin 


1 


X:=0 


2 


for 2 := 1 to A^ do 


3 




for j := —K to K do 


4 




X := X + Un{Zij)Un\zij)l^ij 

end 




end 


5 


X .^ ±,tlsinc^n 


6 


for 2 = 1 to A^ do 


7 




for j = —K to K do 


8 




Un \Zij J '. — Un \Zij) -LUn \Zij J 

end 




end 


end 
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Procedure ConiputeF(u„ (^ij), A„ 



,W 



input : The arrays of values u„ '(zij), k £ 0,d, i <E 1,N, j £ —K, K and A„ , k <E 0,d 



(fc) 



p{d) 



output: The array of values F„ (zij), i e 1, A^, j € —K,K, computed by formula (2.7) 
begin 



for i ;= 1 to A^ do 

for j := —K to K do 

for A; := to d — 1 do 

p('i)|'7 N Tp'^'^Uy \^\, .i/C^Vr ^ 

end 

Fn {zij) ■■= K {zi,j) + ui ^ {zij) {q{zi,j) - q{zi^j)} 



end 



end 



end 
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(0), 



?(d)t 



Procedure ConiputeCorrectionForEigenfunction(u„ (z^ j), u)„(zjj), Fn '{zij)) 



(0), 



Ml 



J/' 



input : The arrays of values u„ (zij), Wn{zij) and Fn {zi^ 

i G MV, J e -K,K 
output: The array of values Un (zij), i E 1,N, j E —K, K 
begin 

Vi := 0; V2 := 0// Initialize auxiliary variables 
for z := 1 to A^ do 

for j := —K to K do 
Xi := 0; X2 := 
for / := -K to K do 



8 

9 

10 

11 
12 



s(-l) 



Ii:=Xi + u\i\zi^i)5)_iF^'^\zi^i)^i^i// Compute integral 



.(^)(c\^U) 



f u);'{0^n'{Od^ by the Stenger's formula (3.14) 



l2:=l2 + Wn{zi^i)6)_iF^'^\zi^i)fii^i// Compute integral 



/ WniOFn'iOd^ tiy the Stenger's formula ( |3.14D 
-1 

end 

Xi := Iihsinc + Vi 

X2 := l2hsinc + V2 

Mn (-2j,j) := Wn{Zi^j)Xi - Un\zij)X2 

end 

Vi := Xi 
V2 := X2 
end 
end 

The described algorithm does not clime to be the most efficient. It is evident that 
the highest accuracy, that can be achieved by increasing the rank r of the FD-method, is 



restricted by the accuracy of the quadrature formulas (3.13), (3.14), used in this algorithm. 



The question about optimal choice of the parameters r and K is still a pressing issue. 
In the next section the numerical results obtained by the Algorithm 1 are stated. In the 



first example we apply quadrature formulas (3.13), (3.14) with K = 500 and in the other 
examples we used K = 350. 
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4 Numerical examples 



Example 1. Let us consider the following Sturm-Liouville problem 



L{u{x)) 



d 
dx 






lim 

a;-s>±l 



(1 - x') 



+ (A -g(x))u(x) = 0, xG(-1,1), 
du{x) 



(4.1) 



dx 







with 



q[x) = X. 



(4.2) 



Using the FD-method, described in [5] and the computer algebra system Maple, it is 
easy to obtain 



A, 



A, 



A, 



A, 



A, 



A, 



A 



0, 

0, 

1 
~6' 

11 

1080' 

0, 

47 
~ 34020' 



0, 



uy{x) 

(2)/ \ 

Uo\x) 
Uo\x) 

(5)/ \ 

Uo\x) 



x/2 



V2: 



X 



4 ' 

V2x^ V2 

24 T2"' 

V2x^ 5V2x 



+ 



288 288 ' 
V2x^ V2x^ 311v^ 
+ 



5760 270 259200' 
V2x^ V2x^ 1181\/2a; 

"172800 ^ 2880 " 518400 ' 
y^x^ \f2x^ \Vl'il\f2x'^ 76967^2 

+ 



7257600 53760 21772800 457228800 
To control the accuracy of the results obtained by the FD-method we use the following 
functionals 






(1 _ a;2) IM^ ^ 



lit \ jYl 



-1 



dx 



(4.3) 



m 



" 1 

1 


' d 


.1 


dx 


-1 


'- 



(1 - x') 



d Un (x) 

dx 



+ \\n -q{x)j M„ (X) 



dx 



(4.4) 



Also we compare the FD-method results with eigenvalues obtained by SLAIGN2 (see. 

0) 

In all examples, presented in this section, we apply the FD-method with uniform mesh 
on the interval [—1,1] and A^ denotes the number of subintervals. The results obtained 
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Table 1: Example 1, the data obtained by SLEIGN2 



n 


^n.sl2 


TOL 


IFLAG 





-0.1576634831? + 00 


0.23950Z? - 


-13 




1 


0.2090760651? + 01 


0.19238Z?- 


-13 




2 


0.6024032351? + 01 


0.20718Z? - 


-08 




3 


0.1201112261? + 02 


0.10472Z?- 


- 12 




4 


0.200064954Z? + 02 


0.70433Z? - 


-13 





Table 2: Example 1, the O-th eigenvalue, FD-method with A^ = 1 



m 


m 

Ao 




nt\x) 




HI 
Vn 




m 

Ao — Ao,s;2 







0.0 


1.0 


0.56 


0.157663483 


1 


0.0 


0.29 


0.22 


0.157663483 


2 


-0.1666666667 


0.24e- 1 


0.36 X 10-1 


0.90031875 X 10-2 


3 


-0.1666666667 


0.18 X 10-1 


0.16 X 10-1 


0.90031875 X 10-2 


4 


-0.1564814815 


0.21 X 10-2 


0.48 X 10-2 


0.11819977 X 10-2 


5 


-0.1564814815 


0.24 X 10-2 


0.23 X 10-2 


0.11819977 X 10-2 


6 


-0.1578630218 


0.30 X 10-^ 


0.78 X 10-3 


0.1995426 X 10-3 


7 


-0.1578630218 


0.41 X 10-3 


0.41 X 10-3 


0.1995426 X 10-3 


8 


-0.1576253633 


0.52 X 10-4 


0.15 X 10-3 


0.381159 X 10-4 


9 


-0.1576253633 


0.79 X 10-4 


0.78 X 10-4 


0.381159 X 10-4 


10 


-0.1576713252 


0.10 X 10-4 


0.31 X 10-4 


0.78460 X 10-5 





Table 3: Example 1, the 0-th 


eigenvalue, FD-method with A^ = 


1 




m 


m 

Ao 




(m 

"o 


{x) 








m 

Ao- 


-Ao,s(2 




51 


-0.15766348313775096746 


0.12 X 


10-16 


0.13 X 10-1*^ 


0.39 


X 10-8 


52 


-0.15766348313775096031 


0.16 X 


10-" 


0.59 X 10-17 


0.39 


X 10-8 


53 


-0.15766348313775096031 


0.33 X 


10-17 


0.33 X 10-17 


0.39 


X 10-8 


54 


-0.15766348313775096218 


0.39 X 


10-18 


0.15 X 10-17 


0.39 


X 10-8 


55 


-0.15766348313775096218 


0.84 X 


10-18 


0.86 X 10-18 


0.39 


X 10-8 


56 


-0.15766348313775096169 


0.11 X 


10-18 


0.40 X 10-18 


0.39 


X 10-8 


57 


-0.15766348313775096169 


0.22 X 


10-18 


0.22 X 10-18 


0.39 


XlO-8 


58 


-0.15766348313775096182 


0.28 X 


10-19 


0.11 X 10-18 


0.39 


XlO-8 


59 


-0.15766348313775096182 


0.58 X 


10-19 


0.61 X 10-19 


0.39 


X 10-8 


60 


-0.15766348313775096178 


0.74 X 


10-2" 


0.29 X 10-19 


0.39 


XlO-8 
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Table 4: Example 1, 


FD-method with N 


= 3 








n 


m 


An 




,(m) 
An 






u^\x) 




7n 




771 

Xn ^Xn,sl2 







15 


-0.1576634831377509617898 


0.3 X 10-21 


0.33 X 10-21 


0.95 X 10-22 


0.39 X 10-s 


1 


15 


2.090760648363956948786 


0.7 X 10-20 


0.13 X 10-20 


0.26 X 10-21 


0.1482 X 10-^ 


2 


15 


6.024031655336352711291 


0.7 X 10-20 


0.94 X 10-21 


0.22 X 10-21 


0.5035 X 10-^ 


3 


13 


12.01112256362987127625 


0.1 X 10-19 


0.24 X 10-20 


0.20 X 10-21 


0.4 X lO-'^ 


4 


11 


20.00649533292656299628 


0.1 X 10-19 


0.31 X 10-19 


0.15 X 10-20 


0.7 X lO-'^ 



10 12 14 





n-n=0 
+ -n=1 
o -n=2 
0-n=3 
• - n=4 



2 



8 10 12 14 



Figure 1: Example 1. The graphs of the functions In I 



,(™) 



(x) 



(left) and Inf ?7„) (right). 
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by the FD-inethod are presented in tables |2]- 14] and figure [T] below. The result obtained 
by the SLEIGN2 are presented in table [Tj 

Analyzing the data in tables [2| [3] and |4| we can conclude that the FD-method with 
A^ = 1 converges much more slower then the FD-method with A^ = 3. This fact is in 



good agreement with the results of Theorem 2A Furthermore, as it follows from table 
|4| the convergence rate of the FD-method increases as the index n of the trial eigenvalue 
increases. Figure [T] illustrates the exponential nature of the FD-method's convergence. 



Example 2. As the second example, let us consider problem (|4.1|) with 

5 



q^x) 



In 



— X 



+ X 



12 I V3 
The results obtained with SLEIGN2 are presented in table [5] below. 



(4.5) 



Table 5: Example 2, the data obtained by SLEIGN2 



n 


\i,sn 


TOL 


IFLAG 





-1.98326983D + 00 


0.46748£» - 08 




1 


0.8551876831? + 00 


0.73426L» - 07 




2 


0.489606686L> + 01 


0.35447L» - 07 




3 


0.10418377015 + 02 


0.40228L» - 07 




4 


0.18816396515 + 02 


0.61329L» - 11 





It is worth to emphasize that the problem under consideration do not satisfy the 
conditions of theorem |2.1 However, from the results presented in table [6] and figure 
[2] it follows that the FD-method have successfully handled this problem as opposed to 
SLEIGN2, which gives results with essential errors (see the leftmost colon in the table |6]). 

Table 6: Example 2, FD-method with N = 24. 



n 


m 


in 




,(m) 
An 










rn 




7n 

An ^An^s/2 







8 


~ 1 .9831442709774408386 


0.99 X 10-" 


0.26 X 10-" 


0.20 X 10-i« 


0.125559 X 10-3 


1 


8 


0.85727032837311800023 


0.47 X 10-15 


0.73 X 10-1*^ 


0.44 X 10-" 


0.20826454 x 10-^ 


2 


8 


4.8939506826799075597 


0.98 X 10-1^ 


0.29 X 10-" 


0.51 X 10-18 


0.2116177 X 10-2 


3 


8 


10.420511296257433545 


0.3 X 10-" 


0.74 X 10-1*^ 


0.13 X 10-16 


0.213430 X 10-2 


4 


7 


18.816396521508987920 


0.11 X 10-16 


0.87 X 10-" 


0.39 X 10-18 


0.2 X 10-^ 



Example 3. Finally, let us consider problem (|4.1|) with 

, . 1 



q[x) 



+ ln 



X + 



X 



(4.6) 



We find out that SLAIGN2 does not handle this problem. But the FD-method does. The 
results obtained by the FD-method are presented in table [7] and figure [3j As before, figure 
[3] confirms the exponential convergence rate of the method. 



20 



8 




-10 

-20 

-30 

-40 



"n^ 


, , , 
m 


, , , 

n- n=Ol 


:X 


V 


+ -n=1 
o - n=2 
0-n=3 
• - n=4 


: ln( iin(x) ) 


■*J 


^ 




Figure 2: Example 2. The graphs of the functions In ( Un (x) J (left) and In ( Vn ) (right) 



10 12 14 16 18 



10 12 14 16 18 





Figure 3: Example 3. The graphs of the functions In I Un (x) 1 (left) and In I Vn) (right) 



Table 7: Example 3, FD-method with N = 12. 



n 


TO 


m 
Xn 




,(m) 
An 






u^\x) 




m 





18 


0.40796999146419634 


0.42 X 10^1^ 


0.20 X 10-15 


0.17 X 10-15 


1 


18 


3.4136861164474333 


0.4 X 10-15 


0.17 X 10-15 


0.16 X 10-15 


2 


14 


6.7759537951814352 


0.2 X 10-15 


0.33 X 10-15 


0.10 X 10-15 


3 


14 


13.323487340142488 


0.2 X 10-14 


0.25 X 10-15 


0.28 X 10-15 


4 


9 


20.8431972121837340 


0.1 X 1015 


0.96 X 10-16 


0.13 X 10-16 
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5 Conclusions 

In the present paper we construct and theoretically justified the generalized algorithm 
of the FD-method for solving the Sturm-Liouville problem for differential equation of 



the second order (1.1), (1.2) with piecewise continuous functional coefficient q{x). As it 
follows from Theorem 2.1 , the generalized FD-method, which uses the piecewise constant 
approximation of the function q{x), can be applied for the approximation of eigenvalues 
and eigenfunctions with any nonnegative index n. The convergence rate of the method 
can be increased by decreasing the value \\q{x) — q{x)\\^ ,^ ^, . In the case when g(x) = 
(this case was considered in 1^) the FD-method (if converges) allows us to calculate the 
approximation to the eigensolution analytically. But in general case, when q{x) ^ 0, the 
analytical calculations are almost always impossible and it is necessary to use numerical 
integration methods, such as sine quadratures and Stenger's formula (see [s], [lo]). 

The problems considered in examples 2 and 3 do not satisfy the conditions of Theorem 
1 1.1 1 because of the unboundedness of the function q{x) on [—1, 1]. However, as it was shown 
in the mentioned examples, the method successfully converges whereas the well known 
in the mathematical world package SLAIGN2 either gives not more then three correct 
numbers after decimal point (example 2) or can not handle the problem at all (example 
3). This examples indicate that the FD-method has a considerable potential which are to 
be investigated in further mathematical works. 
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